/*
 *  DirectDynamic.cpp
 *  EPPI-0.0
 *
 *  Created by Juan Carlos Garcia Orden on 15/11/07.
 *  Modified by Markus Lutz
 *  Modified by Sergio Conde
 *
 */
 
#include "directdynamic.h"

DirectDynamic::DirectDynamic()
:	Analysis(),
	tIni(0.0),
	tEnd(1.0),
	tIncre(1.0),
	linSolver("gauss"),
	nlSolver("newton"),
	maxIter(85),
	tol(1e-10),
	integratorName("bdf-1"), // edit
	beta(0.25),
	gamma(0.5),
	alpha(1.0),
	vTime(), // vtime
	intOptions(), // edit: options for the integrator
	numericTangFlag(0),
	projPenalty(0.0),
	deltaProjTol(1.0e-1)
{	isVelNeeded = true;}

// ***

DirectDynamic::DirectDynamic(ifstream &in, string & na, string &ty, int di, double ti, string &gf, string &gp, double &gs, Simulation *pSim)
:	Analysis(na,ty,di,gf,gp,gs,pSim),
	tIni(ti),
	tEnd(1.0),
	tIncre(1.0),
	linSolver("gauss"),
	nlSolver("newton"),
	maxIter(85),
	tol(1e-10),
	integratorName("bdf-1"), // edit
	beta(0.25),
	gamma(0.5),
	alpha(1.0),
	vTime(),	// edit
	intOptions(), // edit: options for the integrator
	numericTangFlag(0),
	projPenalty(0.0),
	deltaProjTol(1.0e-1)
{
	isVelNeeded = true; 
	string keyword;
	string sysName;
	string loadName;					/** Auxiliary string */
	string loadType;					/** Auxiliary string */
	string bodyName;					/** Auxiliary string */
	string iniName;
	string iniType;
	string iniPointName;
	string line;	
	double value;
	vector<double> iniValues;
	vector<double> iniValuesOmega;
	vector<string> iniValuesSt;
	int i;
	stringstream stnum;
	string pointLabel;
	string iniBodyName;
	string auxs;
	char* dispFile;
	string dispDofs;
	vector<string> vdispPointNames;
	string dispPointsetName;
	string dispLaw;
	Process *procPointer;
	double dispVal;
	string snum;
    vector < string > newtonData; //New options to choose maxIter and tol

	while (in >> keyword && keyword.compare("ENDANALYSIS") )
	{
		if(!keyword.compare("ENDTIME"))
			in >> tEnd;
		else if(!keyword.compare("TIMEINCRE"))
			in >> tIncre;
		else if(!keyword.compare("LINSOLVER"))
			in >> linSolver;
		else if(!keyword.compare("NLSOLVER"))
		{
			getline(in,line);
			istringstream iss(line);
            while (iss >> keyword )
            {
				newtonData.push_back(keyword);
			}
			nlSolver = newtonData[0]; //  As usual the user sets nlSolver, with maxIter and tol being set by default.

			if(newtonData.size() == 2) // The user only sets maxIter or tol and tol or maxIter is set by default
            {
				stnum << newtonData[1];
				stnum >> value;
				if (value>=1)
				{
					maxIter = value;
				}
				else tol = value;				
			}
			if(newtonData.size() == 3) // The user sets maxIter and tol
			{
				stnum << newtonData[1];
				stnum >> value;
				if (value>=1)
				{
					maxIter = value;
					stnum.clear();
					stnum << newtonData[2];
					stnum >> value;
					tol = value;									
				}
				else 
				{
					tol = value;
					stnum.clear();
					stnum << newtonData[2];
					stnum >> value;
					maxIter = value;						
				}				
			}
		}
		else if(!keyword.compare("INTEGRATOR"))
		{
			in >> integratorName;	
			if (!integratorName.compare("newmark"))
			{
				readOptions(in);
				// to be done: idea input file contains additional solver options in the following pattern
				// e.g
				// INTEGRATOR integratorname [ solveroptions = variablename1 value1 variablename2 value2]: 
				// e.g:
				// 
				// INTEGRATOR newmark [ alpha 0.25 	beta 0.5 	STEPSIZEVARIABLE TRUE ]
				// readOptions(); saves options in a map <variablename> <value>; and 
				// is given to solver integrator(options); when the integrator is generated.
				integratorType = 0;	// 0 means standard integrator (newmark, etc.)
			}
			// 1 means midpoint
			if (!integratorName.compare("midpoint"))
			{
				integratorType = 1;
			}
			// 2 means consistent
			if (!integratorName.compare("consistent"))
			{
				integratorType = 2;
			}

			implicitFlag = isImplicit(integratorName);
			if (implicitFlag==false)	lumpedMassFlag=true;
		}
		else if(!keyword.compare("TIMESPAN")) 
			createTimeVec(in);  
		// edit ends here
		else if(!keyword.compare("REMOVELOAD"))
		{
			changeLoads = true;
			/** Read load to be removed.		***** TO BE DONE */
		}
		else if(!keyword.compare("REMOVESYSTEM"))
		{
			haveToRebuild = true;
			/** Read system to be removed.		***** TO BE DONE */
		}
		else if(!keyword.compare("REMOVEBODY"))
		{
			haveToRebuild = true;
			/** Read body to be removed.		***** TO BE DONE */
		}
		else if(!keyword.compare("ADDSYSTEM"))
		{
			haveToRebuild = true;
			/** Read system to be added.		***** TO BE DONE */
		}
		else if(!keyword.compare("ADDBODY"))
		{
			haveToRebuild = true;
			/** Read body to be added.			***** TO BE DONE */
		}
		else if(!keyword.compare("GRAVITY"))
		{
			// DIRECTIONS ARE 1, 2 and 3
			double grav;
			int gravDir;
			in >> grav;
			in >> gravDir;
			// This load has to be sent to all bodies in all processes
			parentSim->setGravity(grav,gravDir);
		}
		else if(!keyword.compare("AUGMENTED"))
		{
			augLagFlag = true;
		}
		else if(!keyword.compare("PROJECTION"))
		{
			projFlag = true;
			in >> auxs;
			if(!auxs.compare("PENALTY"))	in >> projPenalty;
			else if (!auxs.compare("AUTO"))
			{
				// autoProjFlag = 0  No auto projection
				// autoProjFlag = 1  Auto projection based on energy
				// autoProjFlag = 2  Auto projection based on velocity constraint
				in >> auxs;
				if(!auxs.compare("ENER"))
				{
					autoProjFlag = 1;
					in >> auxs;
					if(!auxs.compare("EXACT"))
					{
					    in >> deltaProjTol;
					    autoProjType = "Exact";
					}
					else if(!auxs.compare("EXP"))
					{
					    in >> deltaProjTol;
					    autoProjType = "Exp";
					}
					else if (!auxs.compare("RAT") || !auxs.compare("APPROX") )
					{
					    in >> deltaProjTol;
					    autoProjType = "Rat";
					}
					else if (!auxs.compare("QUAD"))
					{
					    in >> deltaProjTol;
					    autoProjType = "Quad";
					}
					else if (!auxs.compare("LINEAR"))
					{
					    in >> deltaProjTol;
					    autoProjType = "Linear";
					}	
				}
				else if(!auxs.compare("DOTPHI"))
				{
					autoProjFlag = 2;
					in >> auxs;
					if(!auxs.compare("EXACT"))
					{
					    in >> deltaProjTol;
					    autoProjType = "Exact";
					}
					else if(!auxs.compare("EXP"))
					{
					    in >> deltaProjTol;
					    autoProjType = "Exp";
					}
					else if (!auxs.compare("RAT"))
					{
					    in >> deltaProjTol;
					    autoProjType = "Rat";
					}
					else if (!auxs.compare("INV") || !auxs.compare("APPROX") )
                    {
                        in >> deltaProjTol;
                        autoProjType = "Inv";
                    }
				}
			}
		}				
		else if(!keyword.compare("TEMPERATURE"))
		{
			double theta;
			in >> theta;
			// This load has to be sent to all bodies in all processes
			parentSim->setProcTemperature(theta);
		}		
		else if(!keyword.compare("BEGINLOAD"))
		{
			changeLoads = true;		/** Loads change, because a new load is defined next */
			in >> keyword;
			if(!keyword.compare("NAME"))
			{
				in >> loadName;
			} else cout << "*ERROR: keyword NAME is missing or out of place from load \n";
			in >> keyword;
			if(!keyword.compare("TYPE"))
			{
				in >> loadType;
			} else cout << "*ERROR: keyword TYPE is missing or out of place from load \n";
			in >> keyword;
			if(!keyword.compare("BODY"))
			{
				in >> bodyName;
			} else cout << "*ERROR: keyword BODY is missing or out of place from load \n";
			sysName = getStringHead(bodyName);
			procPointer = parentSim->getProcFromSys(sysName);	/** Pointer to the process where the body is */
			procPointer->sendLoadToBody(loadType,sysName,bodyName,in,vTime);	/** Send the load to the body */		
		}
		else if(!keyword.compare("BEGININITIAL"))
		{
			iniValues.clear();
			iniValuesSt.clear();
			
			in >> keyword;
			if(!keyword.compare("NAME"))	in >> iniName;
			in >> keyword;
			if(!keyword.compare("TYPE"))	in >> iniType;
			if(!iniType.compare("position") ||  !iniType.compare("velocity") )
			{
				in >> pointLabel;	// Read the POINT label
				if(!pointLabel.compare("POINT") || !pointLabel.compare("POINTSET"))
				{
					in >> iniPointName;
					sysName = getStringHead(iniPointName);
					getline(in,line);	// I need this to force it to jump to the next line
					getline(in,line);
					//				cout << "line is " << line << "\n";
					istringstream iss(line);
					while (iss >> keyword )
					{
						//					cout << "keyword is " << keyword << endl;
						iniValuesSt.push_back(keyword);
					}
					for(i=1; i<iniValuesSt.size(); ++i)	// Starts from 1 because the first element is the keyword VALUE
					{
						stnum << iniValuesSt[i];
						stnum >> value;
						iniValues.push_back(value);
						stnum.clear();
					}
					procPointer = parentSim->getProcFromSys(sysName);
					if(!iniType.compare("position"))
					{
						if(!pointLabel.compare("POINT"))	procPointer->sendIniPosToPoint(sysName,iniPointName,iniValues);
						if(!pointLabel.compare("POINTSET")) procPointer->sendIniPosToPointSet(sysName,iniPointName,iniValues);					
					}
					if(!iniType.compare("velocity"))
					{
//						cout << sysName << " " << iniPointName << endl;
						if(!pointLabel.compare("POINT"))	procPointer->sendIniVelToPoint(sysName,iniPointName,iniValues);
						if(!pointLabel.compare("POINTSET"))	procPointer->sendIniVelToPointSet(sysName,iniPointName,iniValues);
					}
					iniValues.clear();
					iniValuesSt.clear();
					cout << "Initial condition " << iniName << " set " << endl;
				}
				else if(!pointLabel.compare("BODY"))
				{
					in >> iniBodyName;
					sysName = getStringHead(iniBodyName);
					getline(in,line);	// I need this to force it to jump to the next line
					getline(in,line);
					//				cout << "line is " << line << "\n";
					istringstream iss(line);
					while (iss >> keyword )
					{
						//					cout << "keyword is " << keyword << endl;
						iniValuesSt.push_back(keyword);
					}
					for(i=1; i<iniValuesSt.size(); ++i)	// Starts from 1 because the first element is the keyword VALUE
					{
						stnum << iniValuesSt[i];
						stnum >> value;
						iniValues.push_back(value);
						stnum.clear();
					}
					procPointer = parentSim->getProcFromSys(sysName);
					if(!iniType.compare("position"))
					{
						procPointer->sendIniPosToBody(sysName,iniBodyName,iniValues);
					}
					if(!iniType.compare("velocity"))
					{
//						procPointer->sendIniVelToPoint(sysName,iniPointName,iniValues);
					}
					iniValues.clear();
					iniValuesSt.clear();
					cout << "Initial condition " << iniName << " set " << endl;
				}
			}
			if(!iniType.compare("rigidBodyVelocity") )
			{
				iniValues.clear();
				in >> auxs;		// Read the BODY label
				in >> bodyName;
				sysName = getStringHead(bodyName);
				procPointer = parentSim->getProcFromSys(sysName);

				in >> auxs;		// Read the POINT label
				in >> iniPointName;
				
				getline(in,line);	// I need this to force it to jump to the next line
				getline(in,line);
								//cout << "line is " << line << "\n";
				istringstream iss(line);
				while (iss >> keyword )
				{
								//cout << "keyword is " << keyword << endl;
					iniValuesSt.push_back(keyword);
				}
				for(i=1; i<iniValuesSt.size(); ++i)	// Starts from 1 because the first element is the keyword V0
				{
					stnum << iniValuesSt[i];
					stnum >> value;
					iniValues.push_back(value);
					stnum.clear();
				}
				iniValuesSt.clear();
//				getline(in,line);	// I need this to force it to jump to the next line
				getline(in,line);
				//cout << "line is " << line << "\n";
				istringstream iss2(line);
				while (iss2 >> keyword )
				{
					//cout << "keyword is " << keyword << endl;
					iniValuesSt.push_back(keyword);
				}
				for(i=1; i<iniValuesSt.size(); ++i)	// Starts from 1 because the first element is the keyword OMEGA
				{
					stnum << iniValuesSt[i];
					stnum >> value;
					iniValuesOmega.push_back(value);
					stnum.clear();
				}
				/*
				cout << "iniValues = " << endl;
				printSTLVector(iniValues);
				cout << "iniValuesOmega = " << endl;
				printSTLVector(iniValuesOmega);
				*/
				procPointer->sendIniVelToBody(sysName,bodyName,iniPointName,iniValues,iniValuesOmega);
				iniValues.clear();
				iniValuesOmega.clear();
				iniValuesSt.clear();
				cout << "Initial condition " << iniName << " set " << endl;
			}
			
		}
		else if(!keyword.compare("MASSMATRIX"))
		{
			in >> auxs;
			if (!auxs.compare("consistent"))	lumpedMassFlag = false;
			else
			{
				lumpedMassFlag = true;
//				parentSim->lProcess->createLumpedMass();
			}
		}
		else if(!keyword.compare("BEGINDISP"))
		{
			//vSysNames.clear();
			in >> keyword;
			if(!keyword.compare("NAME"))
			{
				in >> auxs;
			} else cout << "*ERROR: keyword NAME is missing or out of place from load \n";
			in >> keyword;
			if(!keyword.compare("POINTS"))
			{
				getline(in,line);
				//cout << "line is " << line << "\n";
				istringstream iss(line);
				while (iss >> keyword )
				{
					//in >> auxs;
					vdispPointNames.push_back(keyword);
				}
				//cout << "Ha pedido puntos " << endl;
			}
			else if(!keyword.compare("POINTSET"))
			{
				in >> dispPointsetName;
				//cout << "Pide pointset " << endl;
			}
			in >> keyword;
			if(!keyword.compare("DOFS"))
			{
				//cout << "Pide dofs " << endl;
				in >> dispDofs;
			}
			in >> keyword;
			if(!keyword.compare("FILE"))
			{
				in >> auxs;
				dispFile=new char[auxs.size()+1];
				dispFile[auxs.size()]=0;
				memcpy(dispFile,auxs.c_str(),auxs.size());				
				
				//cout << "Ha pedido file " << endl;
				parentSim->createHistory(dispFile);
			}
			else if(!keyword.compare("VALUE"))
			{
				in >> dispVal;
				parentSim->createHistory(dispVal);
				auxs = "value";
				stnum.clear();
				stnum << dispVal;
				stnum >> snum;
				auxs+="-";
				auxs+=snum;
				dispFile=new char[auxs.size()+1];
				dispFile[auxs.size()]=0;
				memcpy(dispFile,auxs.c_str(),auxs.size());				
			}
			// Sending pointset and points
			if(vdispPointNames.size()==0)	// Pointset
			{
				sysName = getStringHead(dispPointsetName);
				procPointer = parentSim->getProcFromSys(sysName);
				procPointer->addImposedDispPointSet(dispPointsetName,dispDofs,dispFile);
			}
			else	// Points
			{
				// Loop on points
				for(i=0; i<vdispPointNames.size(); ++i)
				{
					sysName = getStringHead(vdispPointNames[i]);
					procPointer = parentSim->getProcFromSys(sysName);
					procPointer->addImposedDispPoint(vdispPointNames[i],dispDofs,dispFile);
				}
			}
			vdispPointNames.clear();
		}		
	}

/*
	cout << "Los valores iniciales leidos son " << endl;
	for (i=0; i<iniValues.size(); ++i)	cout << iniValues[i] << endl;
*/	
	if(!nlSolver.compare("newton-numeric")) numericTangFlag=1; else numericTangFlag=0;
//	if(!linSolver.compare("gauss")) lmx::setLinSolverType(0);
//	else if(!linSolver.compare("cg")) lmx::setLinSolverType(2);

}

// ***
/*
DirectDynamic::DirectDynamic(string & na, string & ty, int di, string ls, string nls, string integ, double be, double gam, double al, double tI, double tE, double tIn, bool rebuild, bool reset, string& gf, Process *p)
:	Analysis(na,ty,di,rebuild,reset,gf),
	tIni(tI),
	tEnd(tE),
	tIncre(tIn),
	linSolver(ls),
	nlSolver(nls),
	integratorName(integ),
	beta(be),
	gamma(gam),
	alpha(al)
{ isVelNeeded = true; }
*/
// ***

double DirectDynamic::calculate()
{
	if (parentSim->lProcess==0) cout << "Starting calculation of analysis " << name << ", static ... \n";
	double time = 0; // edit: t added ... maybe better, if it is a local variable of analysis or static...
	int numStep = 0; // edit:
	int i;
	int j;
	int numIter = 0;
	//int maxIter = 30;
	bool itConverged;
// 	double tol = 1.0E-10;
	double minTol = tol*tol;	// Tolerance to account for linear problems, that needs just one iteration
	double resmod ;			// resmod is a relative residual energy.  force residual * q increment between iterations / initial energy
	double resmodMax;
	double resmodIni;
	int numInterfaceIter;
//	numInterfaceIter = (parentSim->nProc)*(parentSim->nProc)+1;	// Number of interface iterations
//	numInterfaceIter = (parentSim->nProc)+15;
	numInterfaceIter = (parentSim->nProc)-1;
	vector<double> deltaq(parentSim->lProcess->numDofs);
	
	// inizialisation of solver
	//double h;	
	//h = vTime[1]-vTime[0]; // cambiar!!! paso de tiempo variable, constante
	time = vTime[0];
	Integrator integrator(integratorName,h,parentSim, intOptions);	// edit: create object solver // 

//	integrator.setPartials(h,inte);
	double partial_dq = integrator.getPartial_dq();
	double partial_ddq = integrator.getPartial_ddq();
	
	// Calculation mode of the residual and tangent, and the convergence check
	// 0: for relative energy residual convergence check
	// 1: for balanced relative residual convergence check
	// 2: for HHT (TO BE DONE)
	int modeCalc = 0;	
	
	parentSim->lProcess->writeToFile("--------------------------");
	parentSim->lProcess->writeToFile("PERFORMING CALCULATION ...");
	parentSim->lProcess->writeToFile("");
	parentSim->lProcess->writeToFile("Number of interface iterations: ", numInterfaceIter);
	parentSim->lProcess->writeToFile("");

	parentSim->lProcess->initializeRes(modeCalc,augLagFlag);
	parentSim->lProcess->initializeTang(modeCalc,projFlag,autoProjFlag, lumpedMassFlag, implicitFlag);

	parentSim->lProcess->createConvergedVectors();
	parentSim->lProcess->setQConv();	// Store converged value of q	
	parentSim->lProcess->update(integratorType);	// Update data at components

	if (implicitFlag==true)	// implicit integrators
	{
		// true means that mass matrix is computed
		parentSim->lProcess->myEvaluationImplicit(time, modeCalc, true); // calculate acceleration in process
	}
	else 
	{
		parentSim->lProcess->myEvaluationExplicit(time, true, lumpedMassFlag); // calculate acceleration in process
	}
	//if (parentSim->lProcess->num==1)	cout << "Acabo de hacer la primera myEvaluationExplicit " << endl;
	
	if(printIterFlag==1 && parentSim->lProcess->num==0)	parentSim->lProcess->createIterOutFile();	//Just for one process, all have the same number of iterations

	if(printStressFlag==1 || plotStressFlag==1)	parentSim->lProcess->fillMSmooth(time);	// Smoothing matrix, just once calculation at the begining

	int timeFlag = 1;

	//***** Roberto
	int plotStep = 0;
	int nShot;

	nShot = (int) (graphicsStep/(vTime[1]-vTime[0]));

	if (nShot < 1 or graphicsStep == 0.0)
	{
			nShot = 1;
	}

	// cout << "nShot= " << nShot << endl;

	for (numStep = 0; numStep < vTime.size()-1; numStep++) 			// loop over all timesteps
	{
        parentSim->lProcess->setQConv();                    // Store converged value of q
		parentSim->lProcess->update(integratorType);        // Update data at components
		if (numStep==0)	parentSim->lProcess->setQIni();
        
		if(printDofsFlag==1) parentSim->lProcess->writeFileDofs(time);			// write dof data for current time
		
		if(printStressFlag==1 || plotStressFlag==1)
		{
			parentSim->lProcess->calcSmoothStress(time);
			if (printStressFlag==1)	parentSim->lProcess->writeFileStress(time);
		}

		if(printEnerFlag==1)
		{
			parentSim->lProcess->calcEner();
			parentSim->collectEner();
			parentSim->writeFileEner(time);
		}

		if(printTempFlag==1)	parentSim->lProcess->writeFileTemp(time);

		if(printEntropyFlag==1)
		{
			parentSim->lProcess->calcEntropy();
			parentSim->collectEntropy();
			parentSim->writeFileEntropy(time);
		}
	
		if(printAMFlag==1)
		{
			parentSim->lProcess->calcAM();
			parentSim->collectAM();
			parentSim->writeFileAM(time);
		}
		if(printConFlag==1)
		{
			parentSim->lProcess->calcCon(time);
			parentSim->collectCon();
			parentSim->writeFileCon(time);
		}
		if(printPProjFlag==1)
		{
			parentSim->writeFilePProj(time, projPenalty);
		}

		if (numStep == plotStep)
		{
			parentSim->lProcess->createGraphOutFile(graphicsFormat,graphicsPath,numStep); 	// create graphic output file for current timestep
			parentSim->lProcess->writeFileGraph(graphicsFormat,graphicsPath,numStep,plotDispFlag,plotProcFlag,plotStressFlag);			// write graphic output file for current  timestep

			plotStep = numStep + nShot;
		}
		//***** Juan Luis 
		// int timeFlag = 1; // is necessary to indicate the implicit integrators that a iteration step has started
                
        if (parentSim->lProcess->num==0) cout << "-------- \n";
        if (parentSim->lProcess->num==0) cout << "--- Timestep " <<  numStep + 1  << " of " << vTime.size()-1;  // print current time step
        if (parentSim->lProcess->num==0) cout << "  (" << vTime[numStep+1] << "s of " << vTime.back() << "s)" <<  "\n"; // print current time step in seconds

        parentSim->lProcess->writeToFile("   ");
        parentSim->lProcess->writeToFile("---- Timestep ", numStep );           // print current time step to file
		
		augIter = 0;
				
		// Lagrange multiplier loop
		do{
			if (augLagFlag==true)
			{
				if (parentSim->lProcess->num==0) cout << "# - Iteration LAMBDA " << augIter << "\n";
				parentSim->lProcess->updateLagMul();	// Update Lagrange multipliers at components
			}
			
			time = vTime[numStep+1]; 

			if (implicitFlag==true) // implicit integrators
			{	
				timeFlag = 1;

				/*
				// Following lines do not improve the accuracy.  Needs more testing
				// Newmark need to evaluate converged acceleration for augmented Lagrangian				
				if (!integratorName.compare("newmark") && augLagFlag==true)	
				{
					// false means that mass matrix is not computed
					parentSim->lProcess->myEvaluationImplicit(time, modeCalc,integratorType,false); // calculate acceleration in process
				}										
				*/

				int state = 0;
				int newstate = 0;
 
				do // loop is necessary for thermo-plastic element
				{

					state = newstate;

					// Newton-Raphson loop
					for(numIter=0; numIter<maxIter; numIter++)
					{		
						itConverged = true;
						if (parentSim->lProcess->num==0) cout << "-------- \n";
						if (parentSim->lProcess->num==0) cout << "----NR Iteration " << numIter + 1 << "\n";
						parentSim->lProcess->writeToFile("----NR Iteration ", numIter+1);
				
						integrator.integrate(time, timeFlag, numStep);

						timeFlag = 0; // only during the first iteration step is set to 1

						// Compute isolated residue
						parentSim->lProcess->formIsolatedResidue(time,numericTangFlag, modeCalc,integratorType,h);
						// Get partial derivative from integrator
						partial_dq = integrator.getPartial_dq();
						partial_ddq = integrator.getPartial_ddq();

						//cout << partial_dq << " " << partial_ddq << "\n"; 
				
						// Compute tangent of isolated processes	
						parentSim->lProcess->formIsolatedTangent(time, partial_dq, partial_ddq, numericTangFlag,integratorType,h);	

						// Reset internal interface contributions
						parentSim->lProcess->resetInterInt();	
				
						// Prepare tangent of internal dofs
						if(numIter==0 && numStep==0)	parentSim->lProcess->factorizeInternalMatrix(0);
						else parentSim->lProcess->factorizeInternalMatrix(1);
					
						for(j=0; j<numInterfaceIter;j++)
							// Loop on the interfaces on each process
						{		
							cout << "		Iteration " << j << " at interfaces in process " << parentSim->lProcess->num << "...\n";
			
							// Set external interface contributions
							// true means that KIB, KBB, etc have to be buit from process before adding new values
							if (j==0)	parentSim->lProcess->setInterExt(true);
							else parentSim->lProcess->setInterExt(false);
							// Set internal interface contributions
							parentSim->lProcess->setInterInt();
						} 
						// Form final residue and tangent of each process
						parentSim->lProcess->formResidue(modeCalc);
						parentSim->lProcess->formTangent();
		
						// Check convergence
						if (parentSim->lProcess->num==0) cout << "	Checking convergence ...\n";

						parentSim->lProcess->writeToFile("");

						// CONVERGENCE CHECK
						if(modeCalc==0)	// Convergence based on total residual energy
						{
	//*						resmod = fabs(((parentSim->lProcess)->res)*deltaq);	// residual energy
							resmod = 0.0;
							for (i=0; i<deltaq.size(); ++i)
							{
								resmod+= ((parentSim->lProcess)->res)[i] * deltaq[i];
							}
							resmod = fabs(resmod);
					
							if (numIter==0)	resmod = 2*tol;						// To force the first iteration
							else	if (resmodIni!=0.0 && resmod>=minTol) resmod/= resmodIni;	// resmod is a relative residual energy		
						}
						if(modeCalc==1)	// Convergence based on balanced residual energy
						{
							resmod = (parentSim->lProcess->eneRes)/(parentSim->lProcess->eneResMax);
							if (numIter==0)	resmod = 2*tol;						// To force the first iteration
						}
						// END OF CONVERGENCE CHECK
				
						MPI::COMM_WORLD.Allreduce(&resmod,&resmodMax,1,MPI::DOUBLE,MPI::MAX);
		
						if (resmod > tol)
						{
							cout << "		Process " << rankk << " has NOT converged, Residue = " << resmod << "\n";
							parentSim->lProcess->writeToFile("		Process ",rankk," has NOT converged, Residue = ",resmod);
							itConverged=false;
						}
						else
						{
							cout << "		Process " << rankk << " has converged, Residue = " << resmod << "\n";
							parentSim->lProcess->writeToFile("		Process ",rankk," has converged, Residue = ",resmod);
							if (resmodMax > tol)		itConverged=false;
						}

						if(itConverged==true) break;
						else
						{
							cout << "	Solving and updating process " << parentSim->lProcess->num << "...\n";
							// Solve linear system and update at each process
							if (numIter==0)	parentSim->lProcess->solveUpdate(deltaq,0,time);
							else parentSim->lProcess->solveUpdate(deltaq,1,time);
					
							if(numIter==0)
							{
	//*							resmodIni = fabs(((parentSim->lProcess)->res)*deltaq);	// Initial residual energy
								resmodIni = 0.0;
								for (i=0; i<deltaq.size(); ++i)
								{
									resmodIni+= ((parentSim->lProcess)->res)[i] * deltaq[i];
								}
								resmodIni = fabs(resmodIni);
							}
						}
					}// End of Newton loop

					newstate = parentSim->lProcess->checkComponents();				
					
					if (newstate > state) // every component agrees with its state?
					{
						// if repetition is necessary the step is set to its initial state
						timeFlag = 1; // 
						parentSim->lProcess->q = parentSim->lProcess->qConv;
						parentSim->lProcess->dq = parentSim->lProcess->dqConv;
					}

				} while(newstate > state); // End of repeat loop

				
				if(numIter<maxIter && itConverged==true)	// Convergence is achieved
				{					
					parentSim->lProcess->writeToFile(" ------------------ \n");
					parentSim->lProcess->writeToFile("** All processes converged at iteration ",numIter+1);
					parentSim->lProcess->writeToFile(" ------------------");
					parentSim->lProcess->writeToFile("Solution process ",rankk);
					parentSim->lProcess->writeToFile(parentSim->lProcess->q);
				}
				else
				{
					cout << "*ERROR: Convergence in one or more processes not achieved at maxIter = " << maxIter << "\n";
				}
				if(printIterFlag==1 && parentSim->lProcess->num==0)	// Just for one process, all have the same iterations
				{
					parentSim->lProcess->writeFileIter(time,numIter+1);
				}			
			}
			else // explicit integrators
			{
				// Compute new position and velocity
				
				integrator.integrate(time, timeFlag, numStep);
				timeFlag = 0; // only during the first time step is set to 1
			
				// Compute acceleration (use myEvaluationExplicit of the parent process)
				parentSim->lProcess->myEvaluationExplicit(time, false, lumpedMassFlag);
			}
			
			if (augLagFlag==true)
			{
				if (augIter==0)	augConvergenceTot=false;
				else 
				{
					augConvergence = true;
					parentSim->lProcess->checkConvLagMul(augConvergence);	// Check convergence of Lagrange multipliers
					MPI::COMM_WORLD.Allreduce(&augConvergence,&augConvergenceTot,1,MPI::BOOL,MPI::LAND);
				}
//				cout << "en directDyn, augCOnvergenceTot es " << augConvergenceTot << endl;
			}			
			augIter+=1; 
		}	while(augLagFlag==true && augConvergenceTot==false && augIter<maxIterLambda);
		// End of Lagrange multiplier loop
		//
		// Mass-matrix projection
		if (projFlag==true)
		{
			if (autoProjFlag!=0)	projPenalty = parentSim->lProcess->calcProjPenalty(deltaProjTol, autoProjType);
									
			// matType=2 (UMF matrices) have problems with re-use of factorized projection matrix
			if (parentSim->lProcess->matType==2) parentSim->lProcess->projectVel(projPenalty, 0);
			else
			{
				if (numIter==0)	parentSim->lProcess->projectVel(projPenalty, 0);
				else parentSim->lProcess->projectVel(projPenalty, 1);
			}	
		}
	}
	//	End of time loop
	parentSim->lProcess->createGraphOutFile(graphicsFormat, graphicsPath, numStep);		// create graphic output file for last timestep
	parentSim->lProcess->writeFileGraph(graphicsFormat,graphicsPath, numStep,plotDispFlag,plotProcFlag,plotStressFlag);			// write graphic output file for last  timestep
	parentSim->lProcess->closeGraphOutFile(graphicsFormat,numStep);	// close output graphic file with specific commands (needed by pvd format)
	if(printDofsFlag==1) 
	{
		parentSim->lProcess->writeFileDofs(time);			// write dof data for current time
		parentSim->lProcess->closeFileDofs();
	}
	if(printEnerFlag==1)
	{
		parentSim->lProcess->calcEner();
		parentSim->collectEner();
		parentSim->writeFileEner(time);
		parentSim->closeFileEner();
	}
	if(printEntropyFlag==1)
	{
		parentSim->lProcess->calcEntropy();
		parentSim->collectEntropy();
		parentSim->writeFileEntropy(time);
		parentSim->closeFileEntropy();
	}	
	if(printTempFlag==1)
	{
		parentSim->lProcess->writeFileTemp(time);
		parentSim->lProcess->closeFileTemp();
	}
	if(printAMFlag==1)
	{
		parentSim->lProcess->calcAM();
		parentSim->collectAM();
		parentSim->writeFileAM(time);
		parentSim->closeFileAM();		
	}	
	if(printConFlag==1)
	{
		parentSim->lProcess->calcCon(time);
		parentSim->collectCon();
		parentSim->writeFileCon(time);
		parentSim->closeFileCon();		
	}		

	if(printPProjFlag==1)
	{
		parentSim->writeFilePProj(time, projPenalty);
		parentSim->closeFilePProj();		
	}						
	
	if(printStressFlag==1 || plotStressFlag==1)
	{
		parentSim->lProcess->calcSmoothStress(time);
		if (printStressFlag==1)
		{
			parentSim->lProcess->writeFileStress(time);
			parentSim->lProcess->closeFileStress();
		}
	}
	
	parentSim->lProcess->writeToFile(" ------------------ \n");
	parentSim->lProcess->writeToFile("Calculation finished \n");

	return tEnd;
}

// ***

bool DirectDynamic::isImplicit(string &name)
{
	if(	!name.compare("ab-1")	||
		!name.compare("ab-2")	||
		!name.compare("ab-3")	||
		!name.compare("ab-4")	||
		!name.compare("ab-5")	||
		!name.compare("cd"  )   ||	
		!name.compare("rk-1")	||
		!name.compare("rk-2")	||
		!name.compare("rk-3")	||
		!name.compare("rk-4"))	return(false);
	else return(true);
}

// ***

// edit starts here
void DirectDynamic::createTimeVec(ifstream &in)
{
		double t_start; 
		double t_end;    
		//double stepsize; 

		vTime.clear(); 
		in >> t_start;		// start time
		//in >> stepsize;		// time increment
		in >> h;		// time increment
		in >> t_end;		// end time

		// cout << t_start << "\t" << stepsize << "\t" << t_end << "\n"; // borrar

/* 
Two possibilities to create the time vector for a given TIMESPAN e.g t = [0.5  1  3]
	1) creates a time vector that does always contain "t_end": vTime = [0.5  1.5  2.5  3]
	2) creates a time vector that uses a constant increment and therefore may not contain "t_end"
		for example here t will create: vTime = [0.5 1.5 2.5 3.5] 
		note that the last element is rounded: e.g [0.5 1 2.9] will create [0.5 1.5 2.5]

	** change the comments to switch between modes	
*/	
	// 1)
		/*
		for (int i = 0; i <= abs(((t_end - t_start)/stepsize)+0.5); i ++) 
		{
			vTime.push_back(t_start+stepsize*i);				// cout << vTime[i] << "\t"; // borrar
		} 
	
		if ((vTime.back()-t_end) > 0) {
			vTime[vTime.size()-1]=t_end;
		}
		else 
		{
			vTime.push_back(t_end);
		}
*/

	// 2)
		/*
		for (int i = 0; i <= fabs(((t_end - t_start)/stepsize)+0.5); i ++) 
		{
			vTime.push_back(t_start+stepsize*i);				 // cout << vTime[i] << "\t"; // borrar
		} 
		 */
		for (int i = 0; i <= fabs(((t_end - t_start)/h)+0.5); i ++) 
		{
			vTime.push_back(t_start+h*i);				 // cout << vTime[i] << "\t"; // borrar
		} 
	
}
// edit ends here


// ***
void DirectDynamic::readOptions(ifstream & in)
{
	string keyword;
	string line;	
	double value;
	
	intOptions["beta"] = 0.25;
	intOptions["gamma"] = 0.5;
	intOptions["rho_inf"] = 1.0;
	
//	cout << keyword << "\n";
	
	getline(in,line);
//	cout << line << "\n";
	istringstream iss(line);
	
	if (iss >> keyword && !keyword.compare("["))
	{
		while (iss >> keyword && keyword.compare("]") )
		{
			if(!keyword.compare("beta"))
			{
				iss >> value;
				intOptions["beta"] = value;
				cout << intOptions.find("beta")->second << "\n";
			}
			else if(!keyword.compare("gamma"))
			{
				iss >> value;
				intOptions["gamma"] = value;
				cout << intOptions.find("gamma")->second<< "\n";
			}
			else if(!keyword.compare("rho_inf"))
			{
				iss >> value;
				intOptions["rho_inf"] = value;
				cout << intOptions.find("rho_inf")->second<< "\n";
			}			
		}
	}
	else {cout << "no options available\n";}
//	cout << keyword << "\n";
	
}
